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SUMMARY 

The numerical performance of a second-order upwind-based TVD scheme and that 
of a uniform second-order ENO scheme for shock capturing are compared. The TVD 
scheme used in this study is a modified version of Liou, using the flux-difference splitting 
(FDS) of Roe and his “superbee” function as the limiter. The construction of the basic 
ENO scheme is based on Harten, Engquist, Osher, and Chakravarthy, and the 2D 
extensions are obtained by using a Strang-type of fractional-step time-splitting method. 
Numerical results presented include both steady and unsteady, ID and 2D calculations. 
All the chosen test problems have exact solutions so that numerical performance can 
be measured by comparing the computed results to them. For ID calculations, the 
standard shock-tube problems of Sod and Lax are chosen. A very strong shock-tube 
problem, with the initial density ratio of 400 to 1 and pressure ratio of 500 to 1, is also 
used to study the behavior of the two schemes. For 2D calculations, the shock wave 
reflection problems are adopted for testing. The cases presented in this report include 
flows with Mach numbers of 2.9, 5.0, and 10.0. 


*Work funded under Space Act Agreement C99066G. 



INTRODUCTION 


For hyperbolic systems of conservation laws, conventional shock capturing schemes 
are known to yield oscillatory solutions near discontinuities. The requirement of mono- 
tonicity has led to the notion of the total variation diminishing 1 (TVD) property that 
provides a mathematical basis for the subsequent development of many TVD schemes. 
One large class of TVD schemes uses flux limiters to control the amount of anti-diffusive 
flux 1-7 . The limiter can be designed so that a conventional non-TVD scheme may be 
modified to satisfy the TVD condition 1 . Formal extensions of these ideas have also 
been made to a variety of problems in multidimensions 8,9 . 

Although “high-order” TVD schemes generally show oscillation-free and crisp 
shock profiles, they degenerate to first-order accuracy at extremum points of the so- 
lution. Also, it was shown that TVD schemes are at most first-order accurate in 
multidimensions 10 . Partly aimed at removing the local restrictions of the TVD schemes, 
Harten, Engquist, Osher, and Chakravarthy recently developed the uniformly accurate 
essentially non-oscillatory (ENO) scheme 11 ’ 12 . Instead of restricting the total variation 
from increasing as in TVD, the ENO scheme permits the variation to possibly increase, 
but only by an amount on the level of local truncation errors. 

The basic ENO-theory contains many desirable properties. The essential feature 
that distinguishes the ENO construction from other shock capturing schemes is the use 
of piecewise polynomials to obtain an essentially non-oscillatory reconstruction of the 
solution from its cell averages. Because the ENO scheme is relatively new, its numerical 
performance has yet to be fully explored. A numerical comparison of the ENO and 
the TVD schemes based on Chakravarthy and Osher 13 is contained in Chakravarthy 
et al 14 , where the numerical experiment was done on a linear scalar equation showing 
the superiority of the ENO over TVD, which clips at the local extremum of the smooth 
solution. Although clearly the ENO scheme is applicable to the approximation of any 
solution, smooth or discontinuous, it seems that the main purpose of its development 
is still to capture shocks. One would like to study further the behavior of the ENO as 
a shock capturing scheme. One of the major tasks is to investigate whether this new 
approach can be extended to produce higher-order multidimensional shock capturing 
schemes. 

In practical applications a uniform second-order scheme is very desirable. In this 
report, we investigate the numerical performance of a second-order upwind-based TVD 
scheme and that of a uniform second-order ENO scheme for shock capturing. We 
understand that not all TVD schemes perform equally well due to the differences in 
construction. The TVD scheme used in this study is a modified version of Liou 7 , using 
Roe’s flux-difference splitting (FDS) 16 and his “superbee” function 5 as the limiter. 
The construction of the basic ENO scheme is based on Harten et al 11 and the 2D 
extensions presented here are obtained by using a Strang-type of fractional-step time- 
splitting method 16 . Results presented in the following sections include both steady 
and unsteady, ID and 2D calculations. For ID calculations, we use the standard 
shock-tube test problems of Sod and Lax. A very strong shock-tube problem, with the 
initial density ratio of 400 to 1 and pressure ratio of 500 to 1, is also used to study the 
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behavior of the two schemes. For 2D calculations, we adopt the shock wave reflection 
problems for testing so that numerical results can be compared with the exact solutions. 
The cases presented here include flows with Mach numbers of 2.9, 5.0, and 10.0. 


CONSTRUCTION OF THE SCHEMES 

In this section we describe briefly the construction of both schemes. 

ENO Scheme: 

Rather than repeating Ref.ll, we describe in the following a step-by-step algorithm 
for implementing a uniform second-order ENO scheme using the reconstruction via 
deconvolution (RD). 

For a single conservation law of the form 


u t + /(«) i = 0, t > 0, 

u(z,0) = u 0 (z) , 

the integration of (1) over [a^_i,z l+ i]x[t n ,t„+i] leads to 


u” +1 = u t n 


At 

Ax 



i 

3 



( 1 ) 

( 2 ) 

( 3 ) 


where 



1 

At 



f [u{x^ L ,t n ~\~ 


( 4 ) 


and tZj* denotes the cell-average of u over [x f _i,x i+ i] at In this formulation, one 
sees that the function u in the integral in (4) has to be reconstructed from the cell- 
averages {uj 1 }. The ENO scheme uses piecewise polynomials to obtain an essentially 
non-oscillatory reconstruction u(x,£ n ), and then uses a local Taylor expansion to obtain 
u(x,f) over [t n ,t n +i]- 

At t rt , assume that we have obtained {v t n } approximating {«"}. The {uj 1 } comes 
from the numerical scheme 


,n+l 


At 

Ax 


(/ i+ * - A-*). 


( 5 ) 


where 

i r At 

/i+i = Xt J 0 f( v ( x *+k' tn + ^)) d/x - 

In order to compute 7,+i. in (6), we first reconstruct v from {vf}. For a second-order 
scheme, the reconstructed v at t n takes the form 

v(x,t n ) = + s,(x - x t ), over each [x t _i,x i+ i], (7) 
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where the slope s, is chosen in the following way. 
For each [a^,x,- + |], if 


i””+i - 2 < + #"-ii < i»r+! - , + »?i, 


l v{Zi + 0) = ^<"" + > 

li v ^ Xi +' “ °) = ^( 3 C, - 4v“ + «?_,); 


else set 


Tt v(x<+0) = 2 

£.(« +1 -b) = 5 i = (, & ,-« i *), 

where + 0) and ^v(xi - o) denote one-sided derivatives at i, from the right 

and the left respectively. Then 

a,=M(^v(x.-°),^„( Il+ 0 )), (8) 

where M denotes the minmod function 


M(p 


’ ?) = {o. 


° ™in(\p\, | 9 |), if Sgn{p) = Sgn{q) = a; 
0, otherwise. 


The above algebraic procedure for the reconstruction of v is illustrated in Fig 1 
In the picture shown, we have |< +1 - 2< + ifcj < K +2 - 2v* +l + v*\. Then we 
construct a quadratic polynomial passing through v t -, and u t+1 . The slopes of the 
two one-sided tangents to this polynomial at x, and x i+l are denoted by £v{ Xi + 0) 

a " d di v ( x *+ 1 -0) respectively. Hence, doing this construction over every [x,*x t - +1 l, we 
will get two one-sided slopes at every point u t . The minmod function M is then used 
to select the final slope. 

Next, we construct the local Taylor expansion Vi(x,t) over \x_i_ x . ilxf£ t , ,1 
using the original PDE (l) and i t (x,t n ) = v(x,t n ) as the initial values. FW a second- 
order scheme, we truncate the second and higher-order terms in the Taylor expansion. 

t*\ , NoW ’ lt ls a,so sufficient a second-order scheme to approximate the integral in 
(6) by using the trapezoidal rule. We obtain 

■^•+2 ~ 2 £ ’ ^«)) "h /( U ( x t + A , tn+l))] . (9) 

Here in (9) each f(v(x i+ i,t)) is then approximated by 

f R {vi{*i+^t),v i+l {x i+k ,t)), 
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where / r (vl,vr) = /(V^Ojvl.vj?)) and F(x/t; t>i,, i/r) denotes the self-similar solution 
to the Riemann problem 


vt + f(v) x = 0, t > 0, 

I \ VL ’ Z< °’ 

l vr , x>0. 

Note that a Riemann solver can be adopted here to obtain the values f R . 
Hence, combining the above we have 


v. n+1 = v? 






( 10 ) 


where 


with 



\ [ f R (»i(*,- + x , *»)»«•+ 1 fo+x , *n)) 
+/ R («f(a: 1 - + x,*n+i).Wi+i(*i+x.*»+i))]. 


> In) — v i + 
Vi+\(x i+ L,tn) = «? +1 
v,(x, + i,t n+ i) = t>" + 

0, + l(x i+ i,<n-fl) = «?+! 


Ax 

2 ’ 

Ax 

Ax 

^“ 5 i+l “ At/ 




(ii) 


and the s,’s come from (8). This scheme is uniformly second-order accurate in the 
pointwise sense. 

For a nonlinear system of conservation laws of the form 


u t + f{ u ) x = 0, t > 0, 
u(z,0) = u 0 (z), 


( 12 ) 

(13) 


where u(x,t) = (uj (x, t ), ..., u m (x, t)) T and / = (/i, we describe briefly the 

procedure using the characteristic reconstruction 11 . Let a t (u), i=l,...,m, denote the 
eigenvalues of the Jacobian matrix 


A{u) 


d A 

du 


so that ai(u) < ••• < a m (u). Let ri (u), r m (u) be the corresponding linearly in- 
dependent right-eigenvectors. Also, let Zi(u), ...,/ m (u) denote the left-eigenvectors of 
A(u) so that li(u)rk(u) = 8^. The k ih characteristic variable is defined to be 

w k = /*u, k — 1, m. 
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Now, suppose we know the approximations {«”} to the cell-averages {u(x t , £ n )}, 

where u(x,, t n ) = (ui(x,-, t n ), u m (x,, t n )) T . For each fixed i and a fixed k, we use 
the following values 

and the procedure described in (7), (8) to reconstruct w k over fx.-_i,x , J at t n . Let 

w (x,t n ) denote this reconstructed polynomial over [x,_ i , x , + 1 ] . Then the vector- 
valued characteristic reconstruction v(x,t n ) will be of the form 

m 

v{x,t n ) = ^2 u;fc (x,tn)rfc(v t n ), over [x t _i,x t+ i], 
k~ 1 

Then we construct the local Taylor expansion «,•(*, t) over [x,-_i,x t - + i]x[t n ,t n+1 ] using 
the original PDE (12) and t/,(x,i n ) = v(x,t n ) as the initial values. 

The resulting scheme for the system (12) will take the same form as (10) and 
(ll), interpreted as vector equations. The values for f ^ are similarly obtained from a 
Riemann solver. For the Euler equations, one can use either an approximate Riemann 
solver (e.g. the one developed by Roe 15 ) or an exact Riemann solver. The results 
presented here are obtained by using an exact Riemann solver outlined in Chorin 17 . 

For 2D computations, we present the results obtained by using the Strang-type 
of fractional-step time-splitting method 16 . The final 2D scheme is formally of second- 
order. Since this method has been well-documented in the literature (see, e.g., Ref. 
8,17), we omit the details. 

Upwind-Based TVD Scheme: 

The construction of the upwind-based TVD schemes consists of basically the fol- 
lowing steps: 

(l) Decomposition by upwinding: The conventional representation of spatial flux 
derivatives, such as second- or higher-order, central or upwind difference, is de- 
composed into parts consisting of a first-order upwind flux difference and the re- 
maining higher-order flux differences, called anti-diffusive terms. For example, one 
can choose the following second-order central difference formula 


fi u )x - 2^ x ^ i+l - f'- 1 ) 

and decompose it into a first-order upwinding term plus an anti-diffusive term as 
follows 

= i (A+r + A_/+ ) + i^j A "( A+ ' + - A+ n, 

where A+ and denote the forward and backward differences, / + and f~ denote 
the positive and negative fluxes. 

(2) Recombination by limiters: Introduce flux-difference limiters, </>+ and <f >~ , to the 
anti-diffusive terms, 

f* = + A “/ + ) + ^ A '(^ +A+ / + - -r a+/“). 

4 
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(3) Time evolution by TVD: Choose a time integration scheme, e.g., the Lax-Wendroff 
or the two-stage Runge-Kutta scheme as used in the present study, and determine 
the functional form of limiters for satisfying the TVD condition 1 . The results 
presented here are obtained by using Roe’s “superbee” function 5 as the limiters. 

Hence the critical ingredients are the upwinding and the choice of limiters. In the scalar 
case, there are no essential variations in defining the argument of the limiter function 
among TVD schemes; the variations are wider in the case of systems of equations. 

In Ref. 7, the upwinding procedure is achieved by the flux-vector splittings (FVS) 
of Steger and Warming 18 , and van Leer 19 . Since the construction of the present TVD 
scheme is done in terms of flux differences representing the spatial derivatives, as shown 
above, the flux-difference splittings (FDS) of Roe 15 and Osher 20 become a natural 
vehicle for upwinding. In fact, replacing the subroutine evaluating the FVS by that for 
FDS was the only change needed in the code. 

The procedure developed for the scalar hyperbolic equation is generalized to the 
one-dimensional system of conservation laws, which can be decoupled into scalar equa- 
tions for each of the characteristic variables in the constant coefficient case. For more 
details on the TVD construction and various forms of the limiter functions for the sys- 
tem, see Ref.7. A formal extension to multidimensional equations is made by treating 
the flux derivative in each direction individually, following the spirit of the directional 
splitting as stated earlier in the construction of the ENO scheme. This treatment, while 
straightforward in implementation, will result in a wider smearing at a discontinuity in 
multidimension. 


NUMERICAL TESTS AND COMPARISON 


We have chosen test problems with known exact solutions to make meaningful 
comparisons. For ID calculations, we show the numerical results for three unsteady 
shock tube problems with shock strengths ranging from moderate to very strong. Here 
ID Euler equations are solved in both schemes. TVD calculations are carried out using 
the flux-difference splitting (FDS) 15 and the “superbee” function 5 as the limiter. In 
Fig. 2, 3, and 4, we show both the computed and the exact solutions after certain 
time steps. For comparison in each figure, the TVD results are shown in the (a)-sets 
of figures and the ENO results in the (b)-sets. The initial conditions are also plotted 
and they are given as follows. 

Fig. 2. Sod’s Problem: 


(1,0,1), x<5 

(0.125,0,0.1), x > 5. 


Fig. 3. Lax’s Problem: 


{p,u,p) = 


(0.445,0.698,3.528), x < 5 
(0.5,0,0.571), x > 5. 
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Fig.4. Strong Shock Problem: 


| (400,0,500), 

l (1,0,1), 


x < 5 
x > 5. 


For the Sod problem, we take Ax = 0.1, At = 0.03, and compare the results after 60 
time-steps. Similarly, for the Lax problem, we take Ax = 0.1, At = 0.017, and 85 
time-steps; and for the strong shock problem Ax = 0.05, At = 0.01, and 90 time-steps. 

For ID computations, the Sod and Lax problems involve only moderate strength 
shocks. Here both the TVD and ENO schemes show very good results. However, 
in comparison with the exact solutions, the TVD results appear to be slightly less 
dissipative than the ENO in describing both the shock and contact discontinuities, 
but less accurate in the head of the expansion fan. The shock strength in the third 
test problem is considerably higher than the first two cases, thereby presenting a more 
severe test on the schemes. Here, the ENO results seem to be better than the TVD 
results, although both schemes badly diffuse the contact discontinuity. In Fig.4, the 
severe rounding of energy-distributions is probably caused by the starting errors in the 
Riemann problem. In Fig.4 (a), the “kink” in the velocity-distribution is caused by the 
flux-differencing in the scheme we used, but can be removed by adding a smoothing 
term. TVD schemes based on initial-value reconstruction such as van Leer’s MUSCL, 
in general, do not have this phenomenon. 

For multidimensional computations, in general, the shock capturing numerical 
schemes are formal extensions of the corresponding ID schemes. In our case, both the 
TVD and ENO schemes are formally extended to 2D problems. We test the extended 
schemes on problems involving a regular reflection of an oblique shock wave from a 
plane wall. Here the 2D Euler equations are solved, with incoming Mach numbers 
ranging from 2.9 to 10.0 and a shock angle of 29 degrees. The physical situation is 
shown in Fig. 5. The TVD scheme is formally applied directly to the steady Euler 
equations, while the ENO scheme is formally applied via the Strang-type fractional 
steps to the unsteady Euler equations until steady solutions are reached. The inflow 
and the interior initial conditions are fully specified with free stream values, and the 
conditions at the top boundary are set to satisfy the shock-jump relations with a 
specified shock angle. The variables at the outflow boundary are extrapolated. At the 
solid wall, the vertical velocity component and the gradient of the other variables are set 
to zero. The computational domains contain equally divided meshes with Ax = 0.067 
and Ay = 0.05. For the test cases with Mach numbers of 2.9, 5.0, and 10.0, the 
computational domains contain respectively 61x21, 74x21, and 83x21 points. 

In the ENO computations, the case with Mach number 2.9 is straightforward. In 
the cases with Mach numbers 5.0 and 10.0, we have to degrade the scheme to the 
first-order Godunov scheme near very strong shocks. The reason for doing this is that 
the application of the same algorithm used in the Mach 2.9 case does not produce 
steady solutions for the higher Mach number cases. What happened was that, after 
the shocks were correctly captured, the computed solution would continue to change 
and convergence could not be reached. Our experiments indicate that the difficulty 
occurs when the Mach number reaches 4.0. We understand that in Glaz, Colella, Glass 
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and Deschambault 21 on a different numerical experiment, a similar reduction of order 
was made to their second-order Godunov scheme in the immediate vicinity of a strong 
shock. 

In our 2D computations, although the results from both extended schemes are quite 
acceptable, they are not as crisp as those ID results. In Fig. 6, we show the pressure 
contours of the Mach 2.9 ENO computation. A similar picture is obtained from the 
TVD computation. Since it is impossible to compare the accuracy of the results using 
these contour pictures, we look at the pressure distributions at the cross sections at 
y = 0 and y = 0.5. In Fig. 7, 8, and 9, we show both the computed and the exact 
steady solutions for the cases with Mach numbers of 2.9, 5.0, and 10.0 respectively. 
In comparison with the exact solution in the Mach 2.9 case, both the TVD and ENO 
results are quite good. ENO captures discontinuities using slightly fewer mesh points, 
but TVD is slightly less oscillatory here. For the cases with Mach numbers 5.0 and 
10.0, our extended ENO results are inferior to the corresponding TVD results. This 
is due to the reduction to the first-order Godunov scheme near very strong shocks. 
Further investigation is undertaken to improve the method of implementation. 
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Fig. 1. Geometric ENO reconstruction 
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Fig. 2. Numerical solutions of Sod’s problem 
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(a) TVD 
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Fig. 2. Numerical solutions of Sod’s problem(continued) 
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Fig. 5. Shock Reflection Problems 
M=2.9, 5.0, 10.0; /?=23.28, 16.99, 13.89 degrees respectively 



Fig. 6. Shock reflection problem, pressure contours from ENO computation 
Mach number-2.9, shock angle=29 degrees 
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Fig. 8. Shock reflection problem, static pressure distribution at y=0 and y=0.5, 
Mach number=5.0, shock angle=29 degrees 


20 



P/P INF H/rtNr 

,-,a un n 3? 6 *f 96 128 160 

L2S 16 u ctt F — i — t-= — i 1 r — i 1 1 


ORIGINAL page is 
OF POOR QUALITY 



o 



(a) TVD 


(b) ENO 


Fig. 9. Shock reflection problem, static pressure distribution at y-0 and y-0.5 
Mach number=10.0, shock angle=29 degrees 
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performance can be measured by comparing the computed results to them. For ID calculations, the standard 
shock-tube problems of Sod and Lax are chosen. A very strong shock-tube problems, with the initial density ratio 
ot 400 to 1 and pressure ratio of 500 to 1, is also used to study the behavior of the two schemes. For 2D 
calculations, the shock wave reflection problems are adopted for testing. The cases presented in this report 
include flows with Mach numbers of 2.9, 5.0 and 10.0. 
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